Tarea 5: Bayesian Inference Part II

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Fecha límite de entrega:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Primera Parte: Preguntas Teóricas
  5. Segunda Parte: Elaboración de Código

Objetivo

a la uuuuultima tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la ultima parte del curso, los cuales se enfocan principalmente en aplicar inferencia bayesiana para generar regresiones lineales y estudiar métodos de obtención de la posterior mas poderosos, como es MCMC. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Videos de las clases:

Documentación:

Primera Parte: Preguntas Teóricas

A continuación, se presentaran diferentes preguntas que abordan las temáticas vistas en clases. Por favor responda cada una de estas preguntas de forma breve, no más de 4 o 5 lineas.

Pregunta 1:

Señale algunos beneficios (no mas de dos) que nos brinda realizar una regresión lineal bayesiana sobre una regresión ajustada por mínimos cuadrados.

En primer lugar, realizar un tratamiento bayesiano de la regresión lineal nos permite incorporar inofrmación a priori, lo que la hace más flexible que la regresión lineal ajustada por mínimos cuadrados. Esto genera que los modelos sean más robustos y menos proclives a hacer overfitting. En segundo lugar, permite realizar una interpretación de la incertidumdre del modelo de una forma mucho más “rica”, en el sentido de que en vez de tener los intervalos de confianza de los coeficientes obtenidos a partir de la sampling distribution, se puede tener una posterior.

Pregunta 2:

A continuación se presenta un modelo de regresión lineal bayesiana:

\[y_i \sim Normal(\mu, \sigma)\] \[\mu = \beta_0 + \beta_1*x\] \[\beta_0 \sim Normal(10,2)\] \[\beta_1 \sim Normal(0,1)\] \[\sigma \sim Uniform(0,50)\]

En base al modelo presentado, responda las siguientes preguntas:

  • Describa que representa cada una de las lineas del modelo.

  • Señale los parámetros que conforman a la regresión lineal.

  • Que objetivo cumple en el modelo lineal tener una distribución \(Normal(0,1)\) en el parámetro \(\beta_1\).

  • Que propiedad de las regresiones lineales nos entrega \(sigma\).

Respuesta:

1 - Descripción del modelo:

\(y_i \sim Normal(\mu, \sigma)\) corresponde a la Likelihood.

\(\mu\) se refiere a la media y en el caso bayesiano no corresponde a un parámetro del modelo, si no que \(\mu\) esta determinado de forma determinista por \(\mu = \beta_0 + \beta_1 \cdot x\). Entonces \(\mu\) se determina una vez que se estiman los parámetros \(\beta_0\) y \(\beta_1\) del modelo. Se aprecia que cada \(\mu\) sigue una relación lineal.

\(\beta_0\) corresponde a un parámetro del modelo que posee un distribución normal con media 10 y desviación estándar 2. \(\beta_0\) es el Intercepto del modelo lineal.

\(\beta_1\) corresponde a un parámetro del modelo que posee una distribución normal con media 0 y desviación estándar 1. \(\beta_1\) es la pendiente del modelo lineal. \(\sigma\) corresponde a la desviación estándar de \(y_i\) y distribuye uniforme(0,50), se puede notar que este sigma es fijo para todas las observaciones lo cual es equivalente al supuesto de homocedasticidad de la regresión lineal frecuentista.

2- Parámetros que conforman la regresión lineal: \(\beta_0\), \(\beta_1\), \(\sigma\)

3 - El objetivo que cumple en el modelo lineal tener una distribución \(Normal(0,1)\) en el parámetro \(\beta_1\). En este caso se utiliza un prior regulzarizador, bien concentrado, Esta regularización se utiliza para prevenir overfitting proveniente de la data asignándole valores extremos a \(\beta_1\).

4 - Cuando se asigna una desviación estándar baja, se le está pidiendo al modelo que nunca asigne pesos demasiado grandes a los betas, previniendo overfitting. Si hay un outlier por ejemplo, el modelo puede ignorar aquellos casos extremos.\(\sigma\) corresponde a la desviación estándar. La propiedad de las regresiones lineales que nos entrega \(sigma\) equivale a la propiedad de homocedasticidad del modelo de regresión lineal estándar.

Pregunta 3:

Explique de forma resumida como funciona el algoritmo de Laplace approximation utilizado para la regresión lineal. Señale el porque existe la necesidad de aplicar este modelo y los supuestos que se realizan con este método.

Respuesta:

El algoritmo de Laplace approximation es una técnica que permite aproximar la posterior con una gaussiana. La idea de ocupar laplace approximation es primero encontrar los valores de máximo a posteriori (encontrar los valores que maximizan el posterior). Esto se puede hacer con el posterior no estandarizado usando algún algoritmo de optimización estilo hill climbing, o sea algoritmos que van recorriendo la función hasta encontrar el óptimo. Luego, se ajusta una distribución gaussiana multivariada, centrada en los valores encontrados. Cuando se usa la Laplace approximation se asume que la posterior conjunta sigue una distribución gaussiana multivariada.

Pregunta 4:

Determine si las siguientes afirmaciones son verdaderas o falsas. Justifique las falsas.

  • [F] El algoritmo de metropolis hasting solo funciona si la probabilidad de saltar de B a A es la misma que de A a B.

  • Se sugiere partir con una misma probabilidad pero luego pueden tener probabilidades diferentes. La probabilidad de saltar de B a A siempre es la misma y se ocupa para estimar otras cosas, en metropolis normal debe ser simétrica, pero en metrópolis hasting no es necesario que esta probabilidad sea simétrica.

  • [F] El algoritmo de Gibbs funciona en cualquier caso.

  • Para una gran cantidad de parámetros (del orden de cien a miles de parámetros) el algoritmo de Gibss no funciona correctamente.

  • [V] El algoritmo de metropolis hasting y de Gibbs son el mismo algoritmo, pero este ultimo es una variante del primero.

Pregunta 5:

El algoritmo de Gibbs es mas eficiente que el de metropolis hasting. ¿Como se logra este efecto? ¿Existe alguna limitante de esta estrategia?

El algoritmo de Gibbs si es mas eficiente que el de metropolis hasting, debido a que utiliza propuestas “más inteligentes”. Este efecto se logra con propuestas que se adaptativas en las que la distribución de los valores de los parámetros que se proponen se ajustan de mejor manera, haciéndolos depender de los valores de los parámetros en el momento. Estas propuestas adaptativas son calculadas por el muestreo de gibbs, y dependen del uso de combinaciones conjugadas entre priors y probabilidades. El algotirmo de Gibbs es altamente aleatorio, pero reduce esta aleatoriedad para ganar eficiencia y explotar el conocimiento de la distribución objetivo. Existen limitantes de esta estrategia, entre ellas se puede mencionar que en algunas ocasiones no vamos a querer utilizar priors conjugados y que es un algoritmo que puede volverse ineficiente al tener demasiados (cientos o miles) de parámetros.

Pregunta 6:

De una ventaja y una desventaja del algoritmo HMC.

Una ventaja del algoritmo de Hamiltonian Monte Carlo es que funciona bien para modelos que posean muchos parámetros (orden de cientos o miles) Una desventaja es que, en comparación los algoritmos Metropolos y Gibbs, es computacionalmente más costoso.

Pregunta 7:

Nombre y explique dos propiedades que son deseables en una cadena de Markov.

Una cadena de markov corresponde a una familia de métodos para hacer muestras de un posterior, sin necesidad de asumir su forma. Una propiedad deseable en una cadena de markov es que se pueda acceder a una función que sea proporcional a la posterior, para poder producir muestras de la posterior. Evaluando el producto de la likelihood y el prior, sin normalizar por la evidencia, MCMC nos permite generar valores random que son representaticos de la distribución de la posterior. Otra propiedad deseable es que el estado actual de la cadena depende solamente de el estado anterior (y no de todos los estados anteriores).

Pregunta 8:

Explique cómo cross-validation, criterios de información y regularización ayudan a mitigar o medir los problemas de underfitting y overfitting.

Cross-validation permite medir como se ajustan los datos al modelo de manera práctica. En particular, dado un dataset y un modelo a entrenar, cross-validation particiona el dataset y entrena el modelo con parte del dataset original y luego evalua el modelo entrenado con la otra parte del dataset. Lo anterior permite generar un mecanismo para medir que tan buena es la capacidad de generalización del modelo dado datos reales (dataset de evaluación) y, así, medir los problemas de underfitting y overfitting.

Criterios de información es una técnica que permite medir los problemas de underfitting y overfitting mediante la creación de métricas (estimadores teóricos) de la out-of-sample deviance. Un ejemplo de métrica es el criterio de información Akali information criteria que es una estimación muy simple de la out-of-sample deviance. En particular, se define la AIC de la siguiente forma: \[ AIC = D_{train} + 2 p\] Donde \(p\) es el número de parámetros del modelo y \(D_{train}\) es a in-sample deviance.

Regularización es un técnica que permite mitigar los problemas de overfiting. En particular, la regularización le otorgar menor peso a parámetros del modelo. Por ejemplo, en una regresión lineal bayesiana con \(\mu = \beta_0 + \beta_1 \cdot x\), definir \(\beta_1 \sim Normal(0,1)\) es un tipo de regularización ya que con dicha distribución el parámetro \(\beta_1\) no se sobre-representa en el modelo y, por tanto, evita el overfitting.

Pregunta 9:

Diseñe una DAG para un problema causal inventado por usted de al menos 5 nodos (puede basarse en el ejemplo de Eugene Charniak) usando Dagitty y considere que la DAG tenga al menos: una chain, un fork y un collider. Muestre con dagitty todas las independencias condicionales de su DAG. Explique las independencias usando las reglas de d-separación.

library(dagitty)
over <- dagitty("dag {
                ser_otaku -> hacer_cosplay;
                ser_otaku -> ver_series;
                ver_series <- tener_acceso_a_internet;
                asistir_a_la_comic_con <- hacer_cosplay
                }")
coordinates(over) <- list(
  x=c("ser_otaku"=2, 
      "hacer_cosplay"=1, 
      "ver_series"=3, 
      "tener_acceso_a_internet"=4,
      "asistir_a_la_comic_con" = 0
      ),
  y=c("ser_otaku"=0, 
      "hacer_cosplay"=1, 
      "ver_series"=1, 
      "tener_acceso_a_internet"=0,
      "asistir_a_la_comic_con" = 3
      ))
plot(over)

Respecot a las independencias condicionales, se tiene lo siguiente:

impliedConditionalIndependencies(over)
## a___ _||_ sr_t | hcr_
## a___ _||_ t___
## a___ _||_ vr_s | sr_t
## a___ _||_ vr_s | hcr_
## hcr_ _||_ t___
## hcr_ _||_ vr_s | sr_t
## sr_t _||_ t___

Más en detalle, el bloque de código anterior muestra las siguiente independecias condicionales: 1 - Asistir a la comic con y ser otaku son condicionalmente independientes dado hacer cosplay. 2 - Asistir a la comic con y ver series son condicionalmente independientes dado ser otaku. 3 - Asistir a la comic con y ver series son condicionalmente independientes dado hacer cosplay. 4 - Hacer cosplay y ver series son condicionalmente independientes dado ser otaku.

1 - Por la regla 1, asistir a la comic y ser otaku están d-conectados pero d-separados dado hacer cosplay.

dconnected(over, "asistir_a_la_comic_con", "ser_otaku")
## [1] TRUE
dseparated(over, "asistir_a_la_comic_con", "ser_otaku", c("hacer_cosplay"))
## [1] TRUE

2 - Por la regla 2, asistir a la comic y ver series están d-conectados pero d-separados dado give hacer cosplay.

dconnected(over, "asistir_a_la_comic_con", "ver_series")
## [1] TRUE
dseparated(over, "asistir_a_la_comic_con", "ver_series", c("ser_otaku"))
## [1] TRUE

3 - Por la regla 1, asitir a la comic con y ver series están d-conectados pero d-separados dado hacer cosplay.

dconnected(over, "asistir_a_la_comic_con", "ver_series")
## [1] TRUE
dseparated(over, "asistir_a_la_comic_con", "ver_series", c("hacer_cosplay"))
## [1] TRUE

4 - Por la regla 1, hacer cosplay y ver series están d-conectados pero d-separados dado ser otaku.

dconnected(over, "hacer_cosplay", "ver_series")
## [1] TRUE
dseparated(over, "hacer_cosplay", "ver_series", c("ser_otaku"))
## [1] TRUE

Segunda Parte: Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.1.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(dplyr)
library(tidyr)

# Para realizar plots
library(scatterplot3d)
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine
# Análisis bayesiano
library("rethinking")
## Loading required package: rstan
## Loading required package: StanHeaders
## rstan (Version 2.21.2, GitRev: 2e1f913d3ca3)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## 
## Attaching package: 'rstan'
## The following object is masked from 'package:tidyr':
## 
##     extract
## Loading required package: parallel
## rethinking (Version 2.01)
## 
## Attaching package: 'rethinking'
## The following object is masked from 'package:purrr':
## 
##     map
## The following object is masked from 'package:stats':
## 
##     rstudent
options(warn = -1)

Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:

install.packages("rethinking")

En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:

install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')

Pregunta 1: Regresión Lineal Bayesiana

El objetivo de esta pregunta es introducirlo en la aplicación de una regresión bayesiana. Con esto, buscaremos entender como calcular una regresión bayesiana en base al “motor” aproximación de Laplace, revisando como se comporta la credibilidad de sus predicciones y como la regresión lineal puede llegar a mutar a aplicar una transformación en el vector \(x\). Para responder esta pregunta centre su desarrollo solo en las clases y las funciones que nos brinda la librería rethinking.

Unos expertos en alometría deciden realizar un estudio de las alturas de unos niños en un colegio, buscando generar un regresor lineal bayesiano capaz de predecir la altura en base al peso de los alumnos. Para realizar este trabajo recopilan los datos table_height.csv, quien posee observaciones fisiológicas de 192 alumnos.

Parte I

En conocimiento los datos recopilados por los expertos, le solicitan realizar la siguiente serie de tareas:

  • Defina un modelo de regresión bayesiana, definiendo sus propios priors. Comente cada una de sus asunciones y señale a través de ecuaciones el modelo. Para definir los priors puede ser interesante la información recopilada en el siguiente link: Priors
  • Ajuste el modelo lineal utilizando el método de Laplace approximation. Estudie a través del comando precis los resultados obtenidos y coméntelos.
  • Gráfique el MAP de regresión lineal obtenida, añadiendo al gráfico el intervalo del \(95\%\) que se tiene sobre la media y los valores predecidos de la altura. Comente los resultados obtenidos y señale si su modelo logra ser un buen predictor de los valores estudiados.

Parte II

En base a los resultados obtenidos, el experto que trabaja con usted le señala que las alturas se suelen modelas con pesos logarítmicos, por lo que le sugiere añadir un logaritmo natural en el vector \(x\) que compone su modelo lineal. Realice nuevamente la regresión utilizando un intervalo del \(95\%\) sobre la media y los valores predecidos de la altura. Comente los resultados obtenidos, señalando si el modelo logra ajustar mejor los valores.

Respuesta:

Inicialmente, se leen los datos asociados al dataset table_height.csv

data_height <- read.csv("table_height.csv")
head(data_height)

Ahora, se define el modelo de la siguiente forma:

\[y_i \sim Normal(\mu, \sigma)\] Se define que los valores \(y_{i}\), es decir, el valor para la \(i\)-ésima altura (la función likelihood del modelo) distribuye de forma normal con media \(\mu\) y desviación estándar \(\sigma\).

Ahora, dado que se quiere realizar una regresión lineal univariada sobre la altura (respecto al peso) el valor de \(\mu\) asociado a \(y_{i}\) se define de la siguiente forma: \[\mu = \beta_0 + \beta_1 \cdot x\]

\(\beta_{0}\) representa el intercepto de la variable \(\mu\) y, dados los datos y la información proporcionada en Priors, se define de la siguiente forma:

\[\beta_0 \sim Normal(60,50)\]

\(\beta_{1}\) representa la pendiente de la variable \(\mu\) y, se define de la siguiente forma para regular el modelo y prevenir que el modelo se ajuste muy drásticamente a los datos:

\[\beta_1 \sim Normal(0,1)\]

Finalmente, \(\sigma\) representa la desviación estándar de todas las alturas y, dados los datos y la información proporcionada en Priors, se define de la siguiente forma: \[\sigma \sim Uniform(0,30)\]

Ahora, dado el modelo definido, se ajusta dicho modelo mediante Laplace approximation.

heigth_reg <- quap(
  alist(
    height ~ dnorm( b0 + b1*weight, sigma ),
    b0 ~ dnorm(60, 50) ,
    b1 ~ dnorm( 0 , 1),
    sigma ~ dunif( 0 , 50)
  ) , data=data_height)

Ahora, para estudiar los resultados de ajustar el modelo se utilzará el comando precis.

precis(heigth_reg, prob=0.95)

En general, el comando muestra los aproximaciones gaussianas para las distribuciones posteriores marginales para cada parámetro del modelo (\(\beta_0\), \(\beta_1\) y \(\sigma\)).

En particular, se tiene que, al \(95\%\), el valor para \(b_{0}\) está entre \(55.73\) y \(61.20\), es decir, el \(95\%\) de la probabilidad posterior asociada a la distribución posterior marginal del parámetro \(\beta_0\) se encuentra entre los valores \(55.73\) y \(61.20\)

Para \(\beta_1\) se tiene que, al \(95\%\), el valor para \(b_{1}\) está entre \(2.57\) y \(2.84\), es decir, el \(95\%\) de la probabilidad posterior asociada a la distribución posterior marginal del parámetro \(\beta_1\) se encuentra entre los valores \(2.57\) y \(2.84\).

Para \(\sigma\) se tiene que, al \(95\%\), el valor para \(\sigma\) está entre \(7.59\) y \(9.28\), es decir, el \(95\%\) de la probabilidad posterior asociada a la distribución posterior marginal del parámetro \(\sigma\) se encuentra entre los valores \(7.59\) y \(9.28\)

Nótese que los resultados anteriores son coherentes con los resultados al ajustar el modelo mediante una regresión lineal no bayesiana. Lo anterior se comprueba mediante la información que proporciona el siguiente bloque de código:

summary(lm(height ~ weight, data = data_height))
## 
## Call:
## lm(formula = height ~ weight, data = data_height)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -31.960  -5.083   1.640   6.048  19.142 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 58.23060    1.40454   41.46   <2e-16 ***
## weight       2.72009    0.06865   39.62   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.481 on 190 degrees of freedom
## Multiple R-squared:  0.892,  Adjusted R-squared:  0.8915 
## F-statistic:  1570 on 1 and 190 DF,  p-value: < 2.2e-16

Para graficar el MAP de la regresión lineal obtenida, añadiendo al gráfico el intervalo del 95% que se tiene sobre la media y los valores predecidos de la altura se tienen que realizar los siguientes pasos:

Se samplean valores para los parámetros del modelos a través de las distribuciones posteriores marginales.

post <- extract.samples(heigth_reg, n=1e4)

Se crea la curva asociado a la regresión utilizando dichos valores sampleados.

weight_seq <- seq(from=3, to=48, by=0.5)
mu_link <- function(weight) post$b0 + post$b1*weight
mu <- sapply(weight_seq, mu_link)
mu_mean <- apply(mu, 2, mean)

Se crea el intervalo del 95% sobre la media.

mu_hdpi <- apply(mu, 2, HPDI, prob=0.95)

Se crea el intervalo del 95% sobre las alturas.

height_weight <- function(weight){
  rnorm(n=nrow(post), mean = post$b0 + post$b1*weight, sd=post$sigma)
}
sim_height <- sapply(weight_seq, height_weight)

height_hpdi <- apply(sim_height, 2, HPDI, prob=0.95)

Finalmente, se grafica el MAP de regresión lineal obtenida, añadiendo al gráfico el intervalo del 95% que se tiene sobre la media y los valores predecidos de la altura.

plot(height~ weight, data=data_height, col=col.alpha(rangi2, 0.9))
lines(weight_seq, mu_mean)
shade(mu_hdpi, weight_seq)
shade(height_hpdi, weight_seq)

Nótese que, según el gráfico, el modelo no se ajusta tan bien porque, para los pesos extremos (es decir, para los pesos menores que 10 y los valores mayores que 30), la curva se aleja mucho de los datos asociados a los a dichos pesos.

Ahora, se realiza el mismo procedimiento anterior pero considerando pesos logarítmicos.

En particular, se crea un nuevo dataset donde se representan los pesos de forma logarítmica.

data_height_2 <- data_height
data_height_2$weight<- log(data_height_2$weight)

El modelo para este caso es el mismo modelo a excepción del parámetro \(\beta_1\) que ahora se define de la siguiente forma: \[\beta_1 \sim Lognormal(0,1)\]

Dado el modelo definido, se ajusta dicho modelo mediante Laplace approximation.

heigth_reg_log <- quap(
  alist(
    height ~ dnorm( b0 + b1*weight, sigma ),
    b0 ~ dnorm(60, 50) ,
    b1 ~ dlnorm( 0 , 1),
    sigma ~ dunif( 0 , 50)
  ) , data=data_height_2)

Para observar los resultados de ajustar el modelo se tiene el siguiente bloque de código:

precis(heigth_reg_log, prob=0.95)

Nótese que, en comparación los resultados obtenidos para el primer dataset, las distribuciones posteriores marginales para los parámetros son distintas.

Del mismo modo, los resultados anteriores son coherentes con los resultados al ajustar el modelo mediante una regresión lineal no bayesiana. Lo anterior se comprueba mediante la información que proporciona el siguiente bloque de código:

summary(lm(height ~ weight, data = data_height_2))
## 
## Call:
## lm(formula = height ~ weight, data = data_height_2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -13.3503  -3.0401   0.2292   2.9634  17.7553 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -32.8566     1.9179  -17.13   <2e-16 ***
## weight       50.5309     0.6757   74.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.679 on 190 degrees of freedom
## Multiple R-squared:  0.9671, Adjusted R-squared:  0.967 
## F-statistic:  5592 on 1 and 190 DF,  p-value: < 2.2e-16

De la misma forma que para el dataset sin pesos logarítmicos, para graficar el MAP de la regresión lineal obtenida, añadiendo al gráfico el intervalo del 95% que se tiene sobre la media y los valores predecidos de la altura se tienen que realizar los siguientes pasos:

Se samplean valores para los parámetros del modelos a través de las distribuciones posteriores marginales.

post <- extract.samples(heigth_reg_log, n=1e4)

Se crea la curva asociado a la regresión utilizando dichos valores sampleados.

weight_seq <- seq(from=0, to=5, by=0.1)
mu_link <- function(weight) post$b0 + post$b1*weight
mu <- sapply(weight_seq, mu_link)
mu_mean <- apply(mu, 2, mean)

Se crea el intervalo del 95% sobre la media

mu_hdpi <- apply(mu, 2, HPDI, prob=0.95)

Se crea el intervalo del 95% sobre los valores predecidos de la altura

height_weight <- function(weight){
  rnorm(n=nrow(post), mean = post$b0 + post$b1*weight, sd=post$sigma)
}
sim_height <- sapply(weight_seq, height_weight)

height_hpdi <- apply(sim_height, 2, HPDI, prob=0.95)

Finalmente, se el MAP de regresión lineal obtenida, añadiendo al gráfico el intervalo del 95% que se tiene sobre la media y los valores predecidos de la altura.

plot(height~ weight, data=data_height_2, col=col.alpha(rangi2, 0.9))
lines(weight_seq, mu_mean)
shade(mu_hdpi, weight_seq)
shade(height_hpdi, weight_seq)

Según el gráfico, el modelo se ajusta mucho mejor que el modelo sin pesos logarítmicos ya que, en general, para cualquier peso, la curva está cerca de los datos, es decir, la curva se ajusta de forma correcta a los datos. Por tanto, la utilización de pesos logarítmicos sí ayuda a ajustar mejor el modelo.

Pregunta 2: MCMC

El objetivo de esta pregunta es lograr samplear, mediante la técnica de MCMC, la distribución gamma.

En general la distribución gamma se denota por \(\Gamma(\alpha,\beta)\) donde \(\alpha\) y \(\beta\) son parámetros positivos, a \(\alpha\) se le suele llamar “shape” y a \(\beta\) rate La densidad no normalizada de una distribución gamma esta dada por:

\[ f(x\mid \alpha,\beta) = \begin{cases} x^{\alpha -1}e^{-\beta x} ~ &\text{ si } x > 0\\ 0 ~&\text{si } x \leq 0 \end{cases} \]

  • Implemente el algoritmo de metropolis hasting utilizando \(q(\theta^{(t)} \mid \theta^{(t-1)}) = \mathcal{N}(\theta^{t-1},1)\), importante su función tiene que poder recibir:
    • La condición inicial \(\theta_{0}\).
    • Cantidad de repeticiones.
    • \(\alpha\) y \(\beta\).
    Escriba el algoritmo sin utilizar implementenaciones de la distribución gamma en r.

De ahora en adelante considere \(\alpha = 5\) y \(\beta = \frac{1}{5}\).

  • Considere \(\theta_{0} = 1\), grafique el histograma que obtiene para distintas cantidad de repeticiones, entre sus pruebas tiene que tener al menos una que tenga del orden de \(10^5\) repeticiones ¿Como cambia la distribución en funcion de las repeticiones?
  • Considere distintos valores de \(\theta_{0}\) y una cantidad de repeticiones grande (del orden de \(10^5\)), grafique las distribuciones que obtenga comente sus resultados ¿es importante la condición inicial del algoritmo?.
  • Utilizando \(\theta_{0}\) y cantidad de repeticiones conveniente (de acuerdo a lo que obtuvo en las partes anteriores) compare con la distribución real. Obs: En esta parte puede utilizar la distribución gamma que tiene implementado r para comparar con lo que usted realizo.

Respuesta:

El siguiente bloque de código contiene las funciones necesarias para implementar el algoritmo de metropolis hasting para samplear la distribución gamma.

fun_gamma <- function(x, alpha, beta){
  if(x > 0){
    return(x^(alpha-1)*exp(-beta*x))
  }
  else{
    return(0)
  }
}
    
metropolis <- function(theta_0, n, alpha, beta){
  values <- rep(0, n)
  current <- theta_0
  
  for (i in 1:n){
    values[i] <- current
    
    proposal <- rnorm(1, current, 1)
    up <- fun_gamma(proposal, alpha, beta)
    down <- fun_gamma(current, alpha, beta)
    r <- ifelse(down==0, 1, up/down)
    
    prob <- min(1,r)
    decision <- rbinom(1,1, prob)
    current <- ifelse(decision==1, proposal, current)
  }
  return(values)
}

Adicionalmente, el siguiente bloque de código contiene la función plot_metropolis que permite graficar el histograma de valores generados al samplear una distribución gamma mediante algoritmo de metropolis hasting.

plot_metropolis <- function(theta_0, n,  alpha, betha){
  metro <- metropolis(theta_0, n, alpha, betha)
  
  ggplot(data=data.frame(x = metro), aes(x=x)) +
    geom_histogram(binwidth = 2) +
    labs(x="Valores generados",
         y="Frecuencia",
         title = paste("Histrograma algoritmo metrópolis para distribución gamma\n", "\nN° iteraciones = ", n, "\nalpha = ", alpha, "\nbetha = ", betha))
    
}

Ahora, tomando \(\alpha = 5\), \(\beta = \frac{1}{5}\) y \(\theta_0 = 1\), se tienen los siguientes histogramas para distintas cantidad de repeticiones.

n_s <- c(1e2, 1e3, 1e4, 1e5, 3e5, 5e5)
for(n in n_s){
  plot(plot_metropolis(1, n, 5, 1/5))
}

Nótese que, mientras mayor es la cantidad de repeticiones, se tiene que el histograma se parece cada más a un histograma con una forma similar a una campana de Gauss (en particular, adelantando los resultados obtenidos, se tiene que, mientras mayor es la cantidad de repeticiones, se tiene que el histograma se parece cada más a un histograma generado a partir de muestrear números aleatorios de una distribución gamma con \(\alpha = 5\), \(\beta = \frac{1}{5}\)).

Respecto a los histogramas que varían el valor inicial \(\theta_0\), se tienen los siguientes resultados:

theta_0_s <- c(1, 2, 3, 5, 10, 20)
for(theta_0 in theta_0_s){
  plot(plot_metropolis(theta_0, 1e5, 5, 1/5))
}

Dado los resultados obtenidos, es claro que el valor inicial \(\theta_0\) no es relevante cuando se tiene un número grande de repeticiones ya que los histogramas convergen a histogramas muy similares (independiente del valor incial \(\theta_0\))

Finalmente, se generan dos histograma:

  • El primer histograma se construirá a partir de los valores generados por el algoritmo metrópolis tomando \(\alpha = 5\), \(\beta = \frac{1}{5}\) y \(\theta_0 = 1\) y un número de repeticiones igual a \(100000\).
  • El segundo histograma se construirá a partir de los valores generados por la función rgamma que permite generar valores aleatorios asociados a una distribución gamma. En particular, para dicha distribución se tomará \(\alpha = 5\), \(\beta = \frac{1}{5}\) y se generarán \(100000\) números aleatorios.
plot_metropolis(1, 1e5, 5, 1/5)

real <- rgamma(1e5, 5, 1/5)
ggplot(data=data.frame(x = real), aes(x=x)) +
  geom_histogram(binwidth = 2) +
  labs(x="Valores generados",
       y="Frecuencia",
       title = paste("Histrograma para distribución gamma\n", "\nN° valores generados= ", 1e5, "\nalpha = ", 5, "\nbetha = ", 1/5))

Nótese los dos histogramas son muy similares y, por tanto, se puede afirmar que el algoritmo desarrollado para samplear sobre la distribución gamma cumple con su objetivo, es decir, dicho algoritmo genera números aleatorios acordes a la distribución gamma para valores dados de \(\alpha\) y \(\beta\).

 

A work by CC6104